home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
C/C++ Users Group Library 1996 July
/
C-C++ Users Group Library July 1996.iso
/
listings
/
v_13_06
/
dwyer
/
randcomb.c
< prev
next >
Wrap
Text File
|
1995-02-10
|
5KB
|
208 lines
/* rand_com[b].c
combination multiplicative random number generator
subtract two random numbers modulo moduli1-1 and shuffle
see
L'Ecuyer, Comm. of the ACM 1988 vol. 31
Numerical Recipes in C, 2nd edition, pp. 278-86
shuffling -- Knuth, vol. II
Copyright (c) 1994, 1995 by Gerald P. Dwyer, Jr. */
#include <time.h>
#include <stdlib.h>
#include <float.h>
#include <assert.h>
#define TESTING
#define TRUE (-1)
#define FALSE 0
void init_rand_comb(long *seed1, long *seed2) ;
long get_init_rand(int) ;
long rand_comb(void) ;
long genr_rand_diff(void) ;
long genr_rand(long a, long x, long modulus, long q, long r) ;
/* first generator */
#define MOD1 2147483563L /* modulus */
#define MULT1 40014L /* multiplier */
/* modulus=multiplier*quotient+remainder */
#define Q1 53668L /* quotient = [modulus/multiplier] */
#define R1 12211L /* remainder */
/* second generator */
#define MOD2 2147483399L
#define MULT2 40692L
#define Q2 52774L
#define R2 3791L
#define MOD_COMB (MOD1-1)
#define MIN_VALUE1 1
#define MAX_VALUE1 (MOD1-1)
#define MIN_VALUE2 1
#define MAX_VALUE2 (MOD2-1)
#define MAX_VALUE ( (MOD1<MOD2) ? MAX_VALUE1 : MAX_VALUE2)
#define EXP_VAL 804307721L
#define GENR1(init_rand) genr_rand(MULT1, init_rand, MOD1, Q1, R1)
#define GENR2(init_rand) genr_rand(MULT2, init_rand, MOD2, Q2, R2)
#define IMPOSSIBLE_RAND (-1)
#define STARTUP_RANDS 16 /* throw away initial draws */
#define DIM_RAND 150 /* size of array shuffled */
static long rand1, rand2 ;
static long rand_num = IMPOSSIBLE_RAND ;
static long rands[DIM_RAND] ;
/* initialize generators with seeds
fill rands[] with initial values */
void init_rand_comb(long *seed1, long *seed2)
{
extern long rand1, rand2 ;
extern long rand_num ;
extern long rands[] ;
int i ;
if (*seed1 <= 0 || *seed1 > MAX_VALUE1)
*seed1 = get_init_rand(MAX_VALUE1) ;
if (*seed2 <= 0 || *seed2 > MAX_VALUE2)
*seed2 = get_init_rand(MAX_VALUE2) ;
/* seed the routine */
rand1 = *seed1 ;
rand2 = *seed2 ;
genr_rand_diff() ;
for (i = 1; i < STARTUP_RANDS; i++) /* throw some away */
genr_rand_diff() ;
/* fill the array of randum number values */
for (i = 0; i < DIM_RAND; i++)
rands[i] = genr_rand_diff() ;
/* initialize rand_num for shuffling */
rand_num = rands[DIM_RAND-1] ;
}
/* get a long initial seed for generator
assumes that rand returns a short integer */
long get_init_rand(int max_value)
{
long seed ;
srand((unsigned int)time(NULL)) ; /* initialize system generator */
do {
seed = ((long)rand())*rand() ;
seed += ((long)rand())*rand() ;
} while (seed > max_value) ;
assert(seed > 0) ;
return seed ;
}
/* generate the difference between random numbers
assumes 0 < rand1 < MOD1
0 < rand2 < MOD2
output 1 <= rand_num <= MOD_COMB */
long genr_rand_diff(void)
{
extern long rand1, rand2 ;
long rand_new ;
rand1 = GENR1(rand1) ;
rand2 = GENR2(rand2) ;
rand_new = rand1 - rand2 ;
if (rand_new <= 0)
rand_new += MOD_COMB ;
assert(rand_new >= 1 && rand_new <= MOD_COMB) ;
return rand_new ;
}
/* generate the next value in sequence from generator
uses approximate factoring
residue = (a * x) mod modulus
= a*x - [(a*x)/modulus]*modulus
where
[(a*x)/modulus] = integer less than or equal to (a*x)/modulus
approximate factoring avoids overflow associated with a*x and
uses equivalence of above with
residue = a * (x - q * k) - r * k + (k-k1) * modulus
where
modulus = a * q + r
q = [modulus/a]
k = [x/q] (= [ax/aq])
k1 = [a*x/m]
assumes
a, m > 0
0 < init_rand < modulus
a * a <= modulus
[a*x/a*q]-[a*x/modulus] <= 1
(for only one addition of modulus below) */
long genr_rand(long a, long x, long modulus, long q, long r)
{
long k, residue ;
k = x / q ;
residue = a * (x - q * k) - r * k ;
if (residue < 0)
residue += modulus ;
assert(residue >= 1 && residue <= modulus-1) ;
return residue ;
}
/* get a random number from rands[] and replace it*/
long rand_comb(void)
{
extern long rand_num ;
extern long rands[] ;
int i ;
/* if not initialized, do it now */
if (rand_num == IMPOSSIBLE_RAND) {
rand_num = 1 ;
init_rand_comb(&rand_num, &rand_num) ;
}
/* rand_num from previous call determines next rand_num from rands[] */
i = (int) (((double)DIM_RAND*rand_num)/(double)(MAX_VALUE)) ;
rand_num = rands[i] ;
/* replace random number used */
rands[i] = genr_rand_diff() ;
return rand_num ;
}
#if defined(TESTING)
/* Test the generator */
#include <stdio.h>
void main(void)
{
long seed1=1, seed2=1 ;
int i ;
init_rand_comb(&seed1, &seed2) ;
printf("Seeds for random number generator are %ld %ld\n",
seed1, seed2) ;
i = STARTUP_RANDS + DIM_RAND ;
do {
rand_comb() ;
i++ ;
} while (i < 9999) ;
printf("On draw 10000, random number should be %ld\n", EXP_VAL) ;
printf("On draw %d, random number is %ld\n", i+1, rand_comb()) ;
}
#endif TESTING